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Abstract 

In this paper we describe a new method of calculation of master integrals based on 
the solution of systems of difference equations in one variable. An explicit example 
is given, and the generalization to arbitrary diagrams is described. As example of 
application of the method, we have calculated the values of master integrals for 
single-scale massive three-loop vacuum diagrams, three-loop self-energy diagrams, 
two-loop vertex diagrams and two-loop box diagrams. 
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In this paper we describe a new method of calculation of master integrals based on 
solution of difference equations. We introduce our method immediately by means of an 
example. Let us consider the integral 

where [d^k] — d^k/n^^'^. By using integration-by-parts[l, 2], one finds the recurrence 
relation 

{x-l)J{x)-{x-l-D/2)J{x-l) = Q . (2) 

This recurrence relation is a first-order linear difference equation with polynomial coeffi- 
cients for the function J. A solution of Eq.(2) can be written in the form of an expansion 
in factorial series 

^hnx-K + s + l) 

_ + 1) / I «i I ^2 \ , . 

r{x-K+l) y"'^ x-K+1^ {x-K + l){x-K + 2)^--- J ' 

where /i, K, and Ug are to be determined. By substituting the expansion (3) in the 
difference equation with the help of the Boole's operators[3] tt and p, which allow one to 
rewrite Eq.(3) as a series of powers of p, one finds the values n — 1 and K — — D/2, and 
the recurrence relation between the coefficients 

sa, = a,_i(s + L>/2)(s + L>/2-l) . (4) 

The first coefficient, oq, can be determined by comparing the large-a; behaviour of the 
series, J{x) ~ fx^x^ao — x'^^^ao, with the large-x behaviour of the integral (1) 



[d^'k] 
(A;2 + 1)^ 



I [d^A;]e-^'=' = , (5) 



from which one infers ao = 1. For large s the coefficients grow approximately as 
s\s^^^, similarly to coefficients of asymptotic expansions in powers of 1/x, but with the 
big difference that the factorial growth is compensated by the gamma function in the 
denominator of Eq.(3), so that the generic term of the sum behaves as g^/^-i-^; a^^^j ^j^g 
series turns out to be convergent ior x > D/2. This is a general feature of expansions in 
factorial series. Furthermore, by using the recurrence relation (2) we can evaluate J{x) for 
values of x outside the domain of convergence of the expansion (3). Let us now consider 
the on-mass-shell integral 

^M-/ (P + i)'g'-2p..) - 

Combining identities obtained by integration-by-parts one finds the recurrence relation 
{x - D)I{x -2) + {2x-D- l)I{x - 1) - 3{x - l)I{x) = (1 - D/2)J{x - 1) , (7) 
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where the function (1) appears in the right-hand side. This is a nonhomogeneous difference 
equation of second order. The general solution of this equation has the form 

I(x) = r]Ji(x) + r]2h{x) + h{x) , (8) 

where Ii and I2 are solutions of the homogeneous equation (obtained by setting the right- 
hand side of Eq.(7) to zero), rji and 772 are arbitrary constants, and I3 is a particular 
solution of the nonhomogeneous equation. As before, we look for solutions in the form of 
expansions in factorial series 

where we set 6g ^ — 6q — 1. By substituting the expansions (9) and (3) in the difference 
equation one finds the values 

III =1, /i2 = - 1/3, =1, 

= - D/2 + 1/2, i^2 = - D/2 + 1/2, isTa = - D/2 , 

and the recurrence relations between the coefficients 6^ 

4s6W -{s- 1/2) (7s -D + l/2)b% + 3{s - 3/2) (s - 1/2)%^^}, = , (10) 

(4s + 2)6f - s(7s -D + 4)6^1 + 3s\s - 1)6^2 = {D/2 - l)(a, - sa^-i) , (11) 

where s = s + D/2 — 1, and = bi = for i < 0. As we will sec immediately, the recurrence 
relation for 6^^-* is not needed. The constants r]i and 772 can be determined by comparing 
the large-x behaviour of the series, /i(x) ~ x^^^'^^^^'^, l2{x) ^ (— l/3)'^a;^'^'''^+^/^ and 
h{x) ~ (1/2 — -D/4)x^'^/^, with the large- a; behaviour of I{x) 

One finds that rji ~ V^/'2 and ^72 = 0, therefore the solution I2 can be discarded. A 
numerical value of rji may be alternatively obtained from Eq.(8) evaluated at x = 0, in 
the form rji = (J(l) - /3(0))//i(0). 

Values of J{x), r]ili{x), h{x) and I{x) calculated for D = 4 — 2e are shown in Table 1, 
with the first two terms of the expansion in e; we also show the number of terms of 
the factorial series needed to obtain the results with 16 exact digits. Values for x < 2 
are obtained by using the recurrence relations (2) and (7). We sec that for small x the 
convergence is slow, so that it is convenient to evaluate the series for a; 3> 1, and to use 
repeatedly the recurrence relations in order to obtain J or I for the given value of x; note, 
however, that the recurrence relation (7) is unstable so that each application increases the 
error on I(x) of a factor |/ii///2| = 3, and therefore a value of x too large is not convenient. 

We mention here that there is an alternative way to represent the functions J, Ii and 
I3 by using the Laplace's transformation [3]: 

1 1 
J{x) = Jdt e-^w{t) , Ij{x) = J dt f'-^Vjit) , J = 1,3 , (13) 
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where w and Vj satisfy the differential equations 

t{l-t)w' + {D/2-t)w = , (14) 

t{t - l)(3i + l)v'j + (3^2 -t{D-l)+2- D)vj = Sj^il - D/2)tw , j = 1, 3 . 

(15) 

Now we want to generalize the method used above to the calculation of master integrals 
of arbitrary diagrams. Let us consider a diagram with Nk loops, internal lines and L 
master integrals. Defined a particular order of the denominators Di, . . . , Dn^, we group 
the master integrals in different sets, such that each set contains the master integrals Bmi 
which have Drr, as first denominator: 



/N 
[d^'h] . . . [d'^kN,]-^^ J, , 0<H<i2<...<in: 



(16) 



where m = 1, . . . , N^, — Nk + I, the indices ii, . . . ,in depend on m and /, and the nu- 
merators N^i in general contain products of powers of some irreducible scalar products 
involving the loop momenta, which cannot be simplified by algebraic identities with one 
of denominators. Integrals Bmi, Bm2, • • • , are ordered by increasing number of denomina- 
tors and, for integrals with the same n, by increasing sum of powers of scalar products 
contained in Nmi- 

Changing the exponent of Dm from 1 to x in each master integral, we define the "master 
functions" 

/AT" 
[d^'h] . . . [d^kr,,] — . (17) 

In general, for each possible choice of Dm equal to one of the denominators, a different 
function I mi is defined; but for x = 1 they all must have the same value Imii^) — Bmi- 
This fact is particularly useful for checking the consistency of the calculations. 

By combining by suitable algorithms [4] the identities obtained by integration- by-parts, 
for each value of m we build a system of linear difference equations in triangular form 

Ri l-l Qik 

J2Pu{x)Iml{x + i) = J2Jl(ljklix)^mk{x+j), / = 1, . . . , , (18) 
i=0 k=l j=0 

satisfied by the functions Imj, with polynomial coefficients pu and qjki, where the right-hand 
side of the l-th equation contains the functions from Imi to Im,i-i- The triangular structure 
is particularly useful for simplifying the numerical solution: the equations are solved one 
at a time, for 1 = 1,2,..., Lm', when the equation for Imi has to be solved, all the functions 
Imi, ■ ■ ■ , Im,i-i in the right-hand side are already known. Therefore each single equation is 
similar to Eq.(2) or Eq.(7), with the difference that, for integrals more comphcated than 
Eq.(6), the order of the equation and the degree of polynomials of coefficients are generally 
larger, and the right-hand side contains the sum of several contributions. The solution of 
each equation of the system (18) has the form 

Imi{x) = Y.r^,lf^{x) + l'^"{x) (19) 
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for integer values of x. Each solution I^^{x) of the homogeneous equation, and the 
particular solution I^^{x) of the nonhomogeneous equation are expanded in factorial 
series 

Values of /Xj, Kj and recurrence relations between coefficients a'^p and a^^^^ can be found 
by using the operator method mentioned in the previous example; the values of /ij^^ and 
come from the already known expansions in factorial series of the nonhomogeneous 
term of Eq.(18). The recurrence relations between a'^J^ and a'^J'^^^ are similar to Eqs.(lO)- 
(11), but in general they are of higher order, and the coefficients of the expansions of 
the functions /^i, ... ,Imi~i appear in the right-hand side of the recurrence relations for 

A comparison between the large-x behaviours of the factorial series (20) and of the inte- 
gral Imi{x) allows one to determine which constants r]j are zero (because the corresponding 
values of jjLj and Kj are not compatible) and which must be calculated. 

Let us now consider the large-x behaviour of the master function Imi{x). Choosing the 
momentum routing such that Djn = kl + ml, introducing hyperspherical polar coordinates 
for the integration over ki and separating the angular and radial part, Eq.(17) becomes 

Imlix) = / l} 'L.. m) , (21) 



r(D/2)J {kl + ml 
where / is given by 

m) = ^ / ino[h) I . (22) 

The path of radial integration i is the straight line [0, oo] for euclidean integrals. For 
non-cuclidean integrals, in order to ensure a correct analytical continuation, the path i 
must be deformed[5, 6], turning around some singularities when the external momenta 
get over some values of threshold in the non-euclidean region; for example, considering a 
one-loop self-energy integral with mi = m2 = 1, the path must be deformed if < —1. 
The integral (21) for large x receives contributions from a neighbourhood of the origin and 
from each singularity of the function /{kf) which is on the path i, turning points included. 

We consider here only the case where the path i is not deformed and all the masses are 
non-zero, so that there is only the contribution from the origin. More general cases will 
be discussed elsewhere. Under this assumption and if, for simplicity, /(O) has a non-zero 
finite value, the large-x behaviour of the master function is given by 

Imiix) ^ {my/'-^x~''/'f{0) . (23) 

Only the solutions I^'^{x) with iij = l/ml and Kj + D/2 equal to a non-negative integer 
have a large-a; behaviour compatible with Eq.(23) (often there is only one such solution). 

/(O) is expressed by an integral over k2, ■ ■ ■ ,k]^^, belonging to a diagram with one 
loop less, obtained by setting /ci = in the original integral; it may be calculated by 
inserting an exponent in a denominator and building and solving new difference equations. 
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If the integral (17) contains on-mass-shell denominators D^^i^ — kf — 2pi ■ ki, ... , 
Dm+iN — kf — 2pN ■ ki, the large- a; behaviour of Imi becomes 



T ('r\ ~ (m'^\D/2-N/2-x -D/2+N/2 



m^f^^g^^ (25) 



if /(O) has a non-zero finite value. Note that the number of dimensions of the integral 
in Eq.(24) is the number of on-mass-shell denominators. An additional relation between 
constants rjj may be found by evaluating Eq.(19) at x = 0, as shown in the example. 

Once the values of the arbitrary constants rjj are obtained, the factorial series (20) must 
be summed. If the coefficients a^p or a^'^^^ grow as s\B^ for large s, with \B\ > 1, the 
corresponding factorial series does not converge for any x. This may occur for particular 
values of masses and external momenta of the master integral. In this case we resort to 
the Laplace's transformation mentioned in Eq.(13), by re-expressing Eq.(20) as 



.,NH 



Ifo^^^ = J^t t^-^yHO^^^ ^ jNH^^^ = Y^ j dt ^(t) , (26) 

J=i 

where v^^ and v^^ satisfy respectively homogeneous and nonhomogeneous differential 
equations obtained by substituting these integral representations in the difference equation; 
the differential equations are conveniently solved by expanding in power series the solutions 
around a number of points selected on the interval of integration. 

Now we consider the application of our method to the calculation of master integrals. 
Considered the number of master integrals and the complexity of the systems of difference 
equations for diagrams with two or more loops, the use of an automated tool is necessary. 
We have written on purpose a comprehensive program, called SYS, which automatically 
determines the master integrals of a diagram, builds the systems of difference equations 
and solves the systems numerically by using expansions in factorial series or integral rep- 
resentations. The program contains a simplified symbolic manipulator, used for solving 
the systems of integration-by-parts identities. Expansion in e of the results is obtained by 
expanding everywhere in e and by manipulating truncated series by means of the arith- 
metical routines implemented in the program; in this way all coefficients of expansions in 
e are obtained in numerical form, divergent terms included. Very high-precision results 
(even 50-100 digits) can be easily obtained, as the calculation of master integrals is reduced 
to the sum of convergent factorial or power series in one variable. 

For illustrating the power of our method, we have calculated the master integrals of the 
diagrams shown in Fig.l for D = 4— 2e in the single-scale case, with all masses equal to one, 
with all the external hnes on-mass-shell and, for box diagrams, the Mandelstam variables 
set to s = t = 1 and u = 2. The number of master integrals is 3, 3, 14, 21, 17, 6, 5, 11 
and 11, respectively for the diagrams (a)-(i) (not counting subdiagrams which factorize, or 
with external lines united). The systems of difference equations between master functions 
generated by the program are formed with 44, 28, 245, 304, 362, 68, 81, 139 and 158 
equations, respectively; at present the program makes no use of the symmetries due to the 
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particular values of masses and momenta in order to simplify or reduce the number of the 
equations. Note that for each subdiagram there may be several different master functions, 
according to which denominator is raised to x. Use of Laplace's transformation is required 
for obtaining the master integrals (e), (g) and (i). As an example, the calculation from 
scratch of all the master integrals of the diagram (c), by using a precision of 38 digits, 
requested about 128 hours of CPU time on a 133 MHz PcntiTim PC; wc stress that at this 
preliminary stage of development we directed our efforts to devise tests and cross checks 
rather than to speed up the program. 

For brevity, we show here only the values of the master integrals with numerator equal 
to one, corresponding to the diagrams shown in figure, leaving the list of the remaining 
master integrals to a further more exhaustive paper [4]. As usual, the results have been 
normalized with the division by Fg = F(l + e) raised to the number of loops of the diagram. 

/(a)F7=^ = 2.404113806319e-^ - 10.03527847977 + 35.94478903214e 

- 119.1503507802e^ + 379.7433345095e^ - 1183.320931551e^ + O(e^) , (27) 

7(b)F7^ = 0.9236318265199 - 1.284921671848e + 2.689507626490e^ 

- 5.338399227511e^ + 10.67136736912e^ + O(e^) , (28) 

7(c)F7^ = 0.2796089232826 - 0.1380294113932e + 0.3194688268113e^ 

- 0.4399664109267e=^ + 0.6650515012166e^ + O(e^) , (29) 

7(d)F7^ = 0.1826272375392 - 0.06746690965803e+ 0.1865462420623e^ 

- 0.2498713405447e^ + 0.3796187113121e^ + O(e^) , (30) 

I{e)V-^ = 0.1480133039584- 0.009263002043238e+ 0.10533085373976^ 

- 0.12242920418466^ + 0.18984804575556^ + 0{e') , (31) 

I{f}Tf = 0.2711563494022 + 0.18339410775146+0.053751010587696^ 

+ 0.014461033684196^ + 0.0007461873722766^ + 0(6^) , (32) 

I{g)V-'^ = 0.173896742268 + 0.18166648769626+0.044408991818326^ 

+ 0.022315473857856^ - 0.0030798104797976^ + 0{e') , (33) 

I{h)Tf = 0.1723367907503 + 0.26795784917116+0.135521127551416^ 

+ 0.044685315328336^ + 0.0084306028274596^ + 0{e^) , (34) 

I{i)T-^ = 0.1036407209893 + 0.21424169329876 + 0.140460686713636^ 

+ 0.044371972363886^ + 0{e^) . (35) 

Of these results, only first two terms of Eq. (27) and the first term of Eq. (28) were already 
known[7, 8]. Remaining terms and other results are new. 
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In conclusion, in this work wc have cstabhshed a connection between recurrence rela- 
tions obtained by intcgration-by-parts and difference equations. Using theory of difference 
equation we have developed a new method of calculation of master integrals which can be 
applied to diagrams with any topology and any number of loops. High-precision results 
expanded at will in e can be obtained easily. We have demonstrated the actual validity 
of the method, as first application, by calculating master integrals in the single-scale case. 
Application of the method to different cases, for example diagrams with different masses, 
external momenta over threshold, or zero masses, which involves the solution of higher- 
order difference equations or more laborious determinations of large- a; behaviours, will be 
described in future papers. 
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Figure Captions 

Figure 1: Diagrams corresponding to the master integrals evaluated. 
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Table 1: Values of J{x), 771/1 (x), hix) and I{x), and number of terms of the series. 
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